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The mixing effectiveness, i.e., the enhancement of molecular diffusion, of a flow can 
be quantified in terms of the suppression of concentration variance of a passive scalar 
sustained by steady sources and sinks. The mixing enhancement defined this way is 
the ratio of the RMS fluctuations of the scalar mixed by molecular diffusion alone to 
the (statistically steady-state) RMS fluctuations of the scalar density in the presence 
of stirring. This measure of the effectiveness of the stirring is naturally related to the 
enhancement factor of the equivalent eddy diffusivity over molecular diffusion, and 
depends on the Peclet number. It was recently noted that the maximum possible 
mixing enhancement at a given Peclet number depends as well on the structure of 
the sources and sinks. That is, the mixing efficiency, the effective diffusivity, or the 
eddy diffusion of a flow generally depends on the sources and sinks of whatever is 
being stirred. Here we present the results of particle-based simulations quantitatively 
confirming the source-sink dependence of the mixing enhancement as a function of 
Peclet number for a model flow. 

Keywords; stirring & mixing; transport processes; stochastic particle dynamics; Brownian 
motion; turbulence 

I. INTRODUCTION 

Mixing by fluid flows is a ubiquitous natural phenomenon that plays a central role in 
many of the applied sciences and engineering. A geophysical example is the mixing of 
aerosols (e.g., CO2 supplied by a volcano, say, or by human activity) in the atmosphere. 
Aerosols are dispersed by molecular diffusion on the smallest scales but are more effectively 
spread globally by atmospheric flows. The density — and density fluctuations — of some 
aerosols influence the albedo of the earth and thus have a significant environmental impact. 



Hence it is important to understand fundamental properties of dispersion, mixing, and the 
reduction of concentration fluctuations by stirring flow fields. 

Various aspects of mixing have been the focus of many review articles [UEIEIIIIEIEIEIIH]- 
At the most basic level, the mixing of a passive scalar can be modeled by an advection- 
diffusion equation for the scalar concentration field with a specified stirring flow field. In 
this work we will focus on problems where fluctuations in the scalar field are generated 
and sustained by temporally steady but spatially inhomogeneous sources. The question of 
interest here is this: for a given source distribution, how well can a specified stirring flow 
mix the scalar field? 

Mixing effectiveness can be measured by the scalar variance over the domain. A well- 
mixed scalar field will have a more uniform density with relatively "small" variance while 
increased fluctuations in the scalar density will be reflected in a "large" variance. We put 
quotes around the quantifiers small and large because the variance is a dimensional quantity 
that needs an appropriate dimensional point of reference from which it is being measured. 

Several years ago Thiffeault et al [9j introduced a notion of "mixing enhancement" for a 
velocity field stirring a steadily sustained scalar by comparing the bulk (space-time) averaged 
density variance with and without advecting flow. Mixing is accomplished by molecular 
diffusion alone in the absence of stirring, which can be quite effective on small scales but is not 
generally so good at breaking up and dispersing large-scale fluctuations quickly. Stirring can 
greatly enhance the transport of the scalar from regions of excess density to depleted regions, 
suppressing the variance far below its diffusion-only value. The magnitude of this variance 
suppression by the stirring — the ratio of the variance without stirring to the variance in 
the presence of stirring — is a dimensionless quantity that provides a sensible gauge of 
the mixing effectiveness of the flow. Different advection fields will have different mixing 
efficiencies stirring scalars supplied by different sources. It is then of obvious interest both 
to determine theoretical limits on mixing enhancements for various source configurations 
and to explore whether those limits may be approached — or even perhaps achieved — -for 
particular flows. 

There have been many studies of stirring and mixing of a scalar with fluctuations 
sustained by spatially inhomogeneous sources and sinks. Some of the earliest are by 
Townsend (TUl [H], who was concerned with the effect of turbulence and molecular diffu- 
sion on a heated filament. He found that the spatial localization of the source enhanced 
the role of molecular diffusivity. Durbin |T2] and Drummond [IS] introduced stochastic par- 
ticle models to turbulence modeling, and these allowed more detailed studies of the effect 
of the source on diffusion. Sawford and Hunt [H] pointed out that small sources lead to 
a dependence of the variance on molecular diffusivity. These models were further refined 
by [151 US] . Chertkov et a/ [HI [Tg [HI |20l [21] and Balkovsky & Fouxon [22] addressed 
the case of a random, statistically-steady source. 

In this paper we study the enhancement of mixing by an advection field using a particle- 
based computational scheme that is easy to implement and applicable to a variety of source 
distributions. The idea is to develop a method that accurately simulates advection and 
diffusion of large numbers of particles supplied by a steady source, and to measure density 
fluctuations by "binning" the particles to produce an approximation of the hydrodynamic 
concentration fleld. Unlike a numerical PDE code, a particle code does not prefer speciflc 
forms of the flow or the source (PDE methods generally work best with very smooth fields). 
There is, however, no free lunch: the accuracy of the particle code is ultimately limited 
by the finite number of particles that can be tracked. The limitation to finite numbers 
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of particles inevitably introduces statistical errors due to discrete fluctuations in the local 
density and systematic errors in the variance measurements due to binning. But these 
problems are tractable, and as we will show, the method proves to be quantitatively accurate 
and computationally eflicient for some applications. 



II. THEORETICAL BACKGROUND 

In this section we review basic facts about the mixing enhancement problem as formulated 
by Thiffeault, Doering & Gibbon et al [9J and developed by Plasting & Young [23], Doering 
& Thiffeault j24j, Shaw et al [25j, and Thiffeault and Pavliotis [26j. The dynamics is given 
by the advection-diffusion equation for the concentration of a passive scalar p(t, x) with 
time-independent but spatially inhomogeneous source fleld S{x): 

^ + u-Wp = KAp + Six), (1) 

where k is the molecular diffusivity and u{t, x) is a specifled advection fleld that satisfles 
(at each instant of time) the incompressibility condition 

V ■ u = 0. (2) 

For simplicity, the domain is the (i-torus, i.e., [0, with periodic boundary conditions. We 
limit attention to stirring flelds that satisfy the properties of statistical homogeneity and 
isotropy in space deflned by 
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Ui{-,x) = 0, Ui{-,x)uj{-,x) = —6ij (3) 

where the overbar denotes time-averaging and U is the root mean square speed of the velocity 
fleld, a natural indicator of the intensity of the stirring. These are statistical properties of 
homogeneous isotropic turbulence on the torus, but they are also shared by many other 
kinds of flows. 

We are interested in fluctuations in the concentration p so the spatially averaged back- 
ground density is irrelevant. It is easy to see from ([T]) that the spatial average of p grows 
linearly with time at the rate given by the spatial average of S. Hence we change variables 
to spatially mean-zero quantities 

9{t, x) = p{t, d'^' /^(^' (4) 



and 

six) 



S{x) - -1 1 dV S{x') (5) 



that satisfy 

f)f) 

+ = kAO + s{x). (6) 

(We must also supply initial conditions for p and/or 6 but they play no role in the long-time 
steady statistics that we are interested in.) 
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The "mixedness" of the scalar may be characterized by, among other quantities, the 
long-time averaged variance of p, proportional to the long-time averaged norm of 9, 

{6') := ^lim f£^t^,J d'a: O^t, x) (7) 

The smaller (6'^) is, the more uniform the distribution. The "mixing enhancement" of a 
stirring field is naturally measured by comparing the scalar variance to the variance with 
the same source but in the absence of stirring. To be precise, we compare {O"^) to {0^) where 
9q is the solution to 

r)f) 

= kAOo + six) (8) 

(with, say, the same initial data although these will not affect the long-time averaged fluc- 
tuations). The dimensionless mixing enhancement factor is then defined as 



This quantity carries the subscript because we can also define multiscale mixing enhance- 
ments 125] by weighting large/small wavenumber components of the scalar fluctuations: 



p = (1") 

As discussed in Doering & Thiffeault [21], Shaw et al [25] and Shaw [27], £±i provide a gauge 
of the mixing enhancement of the flow as measured by scalar fluctuations on relatively small 
and large length scales, respectively. We refer to them as enhancement factors because if 
one were to define an effective, eddy, or equivalent diffusivity Ke,p as the value of a molecular 
diffusion necessary to produce the same value of (|V^^^p) with stirring, then Kf.^p = nSp. In 
this paper, however, we will focus exclusively on Sq, the mixing enhancement at "moderate" 
length scales. 

There is a theoretical upper bound on Sq valid for any statistically stationary homoge- 
neous and isotropic stirring field [211 EHl EZ] : 



E,^o\m\vk^ ^^^^ 



where s{k) are the Fourier coefficients of the source and the Peclet number, 

Pe := UL/k, (12) 
is a dimensionless measure of the intensity of the stirring. Generally, we anticipate that 



So is an increasing function of Pe and the estimate in (11) guarantees that £o(Pe) < Pe as 
Pe — > cxD, the "classical" scaling necessary if there is to be any residual variance suppression 
in the singular vanishing diffusion limit. That is, if ^o(P6) ~ Pe then Kg.o has a nonzero limit 
as ft: — *• with all other parameters held fixed. It is natural to refer to any of the possible 
sub-classical scalings as "anomalous" . 
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The upper limit to the mixing enhancement in (11) depends on the stirring field only 
through U via Pe, but it depends on all the details of the source distribution. As studied 
in depth in references [211 12S1 [2Z] , the structure of the scalar source can have a profound 
effect on the high Pe scaling of Sq, notably for sources with small scales. It is physically 
meaningful to consider measure-valued source-sink distributions, like delta-functions, with 
arbitrarily small scales. It is precisely this source size dependence of So(Pe) that motivates 
the development of a computational method that can handle singular source distributions. 

In this study, for computational simplicity and efficiency, we utilize the "random sine 
flow" as the stirring field. In the two-dimensional case this is defined for all time by 

\wsm{27iy/L + 4))i, nT<t<nT + \T] 

u(t,x) = < , , „ 1 , , (13) 

^ ^ \wsm{2nx/L + (j)') j, nT + \T < t < {n + l)T, ^ ' 

where T is the period, n = 0, 1, 2, . . . , and and 0' are random phases chosen independently 
and uniformly on [0, 27r) in each half cycle, which assures the homogeneity of the fiow field. 
In this case, w = \/2U. In the three-dimensional case, we employ 

{w [sin a sin {2'ny/L + 02) + cos a sin {2Ttz/ L + ^3)] i , nT < t < nT + \T] 
w [sin a sin (27tz / L + ^3) + cos a sin (27rx/L + 0i)]j, nT +lT <t <nT +lT; 
w [sin a sin {2ttx/L + ^pi) + cos a sin {2Tcy/L + ^/'2)] k, nT + ^T < t < {n + 1)T, 

where again w = V2U, n = 0,1,2, ... , and a, 0i,2,3 and V'1,2,3 are uniform random numbers 
in [0, 27r) chosen independently every T/3. The angle a randomizes the shear direction to 
guarantee isotropy of the fiow. 



III. NUMERICAL METHOD 

In a particle code for solving the advection-diffusion equation, the concentration field p 
is represented by a distribution of particles. Particles are introduced by generating random 
locations using the properly normalized source S{x) as a probability distribution function, 
then they are transported by advection and diffusion. The particle density, p{t,x), is mea- 
sured by covering the domain with bins counting the number of particles per bin. 

A discrete particle method is employed because it can easily deal with small-scale sources 
such as 5 functions. It is also straightforward to implement with any advection field. The 
downside of a particle method is that it necessarily involves two kinds of errors: the number 
density of particles calculated by dividing the domain into bins is only resolved down to the 
size of the bins, and the measurement of p always includes statistical errors due to the use 
of finite numbers of particles. 



A. Time evolution 



At each time step the system is evolved by advection, diffusion, the source, and sinks. 
An advection-only equation would be solved by moving particles along characteristics, and 
a diffusion-only equation would be solved by adding independent Gaussian noises to each 
coordinate of each particle. With both advection and diffusion we need to solve a stochastic 
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differential equation to determine the proper displacement of the particles during a time 
step. The stochastic differential equation is 



dX = u{t, X)dt + dW 



(14) 



where W{t) is a standard vector- valued Wiener process. 

In order to solve (14), we will consider cases where the displacement due to the noise in 



a subinterval of length T /d (where d is the dimension) is much smaller than the wavelength 
of the random sine flow. This condition is realized better and better as Pe increases. Then, 
during each subinterval, the drift field u{t,X) experienced by each particle can be approx- 
imated by a steady flow with a linear shear. In 2D, for the first half of the period for a 



particle starting at (xo,i/o) = {X{t = 0),y(t = 0)) we approximate (14) by 



2ti 



(Y -ya)dt + ^j2KdWi, 
L (15) 



dX = w sin (27r?/o /L + (j)) dt + w cos (27r?/o /L + 
dY = V2KdW2, 

and for the second half of the period, starting from (xq, i/'q) = {X(t = T/2), Y{t = T/2)), 
dX = V2^dWi, 

dY = w sin {2ttx'q/ L + 0') dt + w cos {2t:x'q/ L + (j)') —(X — x'Q)dt + v2KdW2 ■ 

ij 

Therefore, during the first half period we evolve the position of a particle through a time 
interval At (where At < T/2 need not be small) by the map 



Xq 

yo 



xq + w sin {2nyo/L + (p) At + Ri, 
yo + R2, 



(17) 



where Ri and R2 satisfy 



dRi = S2R2dt + V2^dWi {S2 := 2mvL~^ cos {2'Kyo/L + 0)) 
dR2 = V2KdW2. 



lS^^Kt^ + 2Kt S2Kt^ 



(18) 



The variance-covariance matrix of Ri and R2 is 

EiR;") E{R^R2 
E{R2Ri) E{R^^) 



SoKt^ 



2Kt 



which is realized by 



Ri 
R2 



\S2^Kt^ + 2Kt X A^i + ^J\S2^Kt^ X N2, 
2Kt X A^2, 



(19) 



(20) 
(21) 



where Ni and N2 are independent A^(0, 1) random variables (normally distributed with 
mean and standard deviation 1). The matrix (19) describes the evolution of a passive 
scalar field in a shear flow [2HI EHl ED] . 
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(a) 



(b) 



FIG. 1: (a) A circular Gaussian distribution of particles transported and sheared into (b) an 
elliptical Gaussian cloud. 

Therefore the time evolution map during the first half period is 



xo Xo + wsm{27ryo/L + (l)) At + ^ IS^'^niAt)^ + 2KAt + ^ ^S^^niAt)^ N2 , (22a) 
yo ^ yo + V2KAt N2 . (22b) 

A similar map is employed during the second half of the period. These stochastic maps 
include the shear — in the approximation that the shear remains constant for each particle 
during each half cycle — that causes a "distortion" of a Gaussian cloud of particles; see 
Fig.[Tj 

The same calculations apply for the three-dimensional case: for the first subinterval, the 
time-evolution map is 

Xq Xq + w [sin a sin {27iyQ/L + ^2) + cos a sin {2'nzQ/L + (p^)] 

+ ^J\{S^^ + S.^^)K{Atf + 2KAtN^ + S2^l\K{At)^ N2 + 53y|^(At)i N^, (23a) 

yo^y^ + V2^tN2, (23b) 
zq + V2KAt Ns, (23c) 

where 5*2 = 2itwL^^ smacos{2TTyo/L + (j)2), S3 = 2t:wL~^ cos a cos{2ti zq / L + (p^) , and Ni,N2 
and As are independent A^(0, 1) random variables. The maps for the other subintervals can 
be obtained by cyclic permutation of the coordinates. 

The steady scalar source is realized by introducing a new particle one by one using 
normalized S{x) as a probability distribution function. Numerically, such a probability 
distribution function can be realized by mapping uniform random numbers over [0, 1] with 
the inverse of the cumulative probability distribution function in question. 

New particles are added constantly so the total number of particles continues to increase, 
which slows down the computation. To cope with increasing particles, we implement a 
particle subtraction scheme. Particles eventually get well mixed and "older" particles do 
not contribute to the value of the hydrodynamic variance. There is no added value in keeping 
track of particles that have been in the mix for a very long time, and we can simply remove 



them from the system after a sufficiently long time. It is very important to keep track of 
the "age" of each particle, however, and to only remove sufficiently old well-mixed particles. 
(For example, if a random fraction of particles is removed at regular time intervals, then the 
simulation becomes one of a system particles with a random finite lifetime, described by an 
advection-diffusion equation with an additional density decay term.) 

In order to determine how old particles must be in order to safely remove them without 
affecting the hydrodynamic variance, prior to a full simulation run a test is performed as 
follows. Starting from an initial set of Ni particles located in space according to the source 
distribution, the flow and diffusion are allowed to act and the variance of the number of 
particles per bin, which decays with time, is monitored. The number is of the order of 
the number of particles that are introduced in the full simulation during, say, an interval 
of length T characteristic of the random sine flow. The variance does not decay all the 
way to zero, however, but rather to the variance expected when Ni particles are randomly 
distributed among the bins. The time when the variance achieves this random-distribution 
variance, measured beforehand for a given flow and diffusion strength, is then the required 
"aging" time before particles can be safely removed in the full simulation with the steady 
source. Such a trial run is performed for each flow, diffusion strength, source distribution 
and particle number because this "mixing time" depends on all these factors. Futher details 
of the criteria for removing old particles and extensive tests and benchmark trials may be 
found in Ref. [HT] . 



B. Variance calculation and background noise 

The variance {6"^) is measured by monitoring the fluctuations in the number of particles 
per bin, and time-averaging. In d dimensions the domain is divided into l'^ bins and the 
code calculates {n"^), where n is the number of particles in a bin. Then {6"^) is initially 
approximated by 

{n')-{n)' = {L/lf'{9'). (24) 

We say "initially" because the expression above includes both the hydrodynamic fluctuations 
of interest and discreteness fluctuations resulting solely from the fact that each bin contains 
a flnite number of particles. 

The subtraction scheme eliminates the "well-mixed" particles that do not contribute to 
the value of the hydrodynamic variance. But even if the system were completely mixed so 
that theoretically, {6"^) = 0, the measured variance (n^) — (n)^ would be (very close to, for 
small bins) (n), which is on the order of N/l'^, where is the total number of particles in 
the domain. This follows from the fact that 6{t, x) is represented in this particle method by 



only a flnite number of particles in each flnite size bin. That is, {6'^) as deflned by (24) is 
nonzero even when the particles are uniformly distributed: then the bulk variance includes 
fluctuations as if N particles were randomly thrown in l'^ bins. The helpful fact is that 
the bulk variance contribution from these background fluctuations due to flnite numbers of 
particles in the bins does not depend on (i.e., is uncorrelated with) the hydrodynamic density 
variation from bin to bin. The total contribution to the variance is the sum of the "extra" 
variance in each bin which is linear in the (mean) number of particles in each bin. Hence 
the sum of the variances is ~ and the bulk variance contribution from the background 
fluctuations, Nl'^/L'^'^, can simply be subtracted from the initial estimate for {6"^) in (24). 
The net result is our measured value of the hydrodynamic variance. 
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FIG. 2: Square-shaped source of two different 
source distribution over the squares. 
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In addition to the inevitable fluctuations due to discreteness, density variations are ob- 
served only down to the length scales ~ L// because of the binning density, which is another 
source of error in this procedure. We use I > 100, which tests and benchmark studies 
indicate is sufficient for the examples studies here pij . 

The variance is calculated once for each subinterval, and the instant when it is calculated 
is determined randomly in order to obtain an unbiased time average. Thus each subinterval 
is divided into two parts, before and after variance calculation, and the particle transport 
and source processes are appropriately adapted. The final measured quantities are long time 
averages that are observed to be converged to within the error indicated on the plots below. 



IV. RESULTS 



In order to investigate the effect of source-sink scales on maximal and actual mixing 
enhancements, we performed a series of simulations for square-shaped sources of various 
sizes a < L as illustrated in Fig. |2j 

Fig 3(a) shows the upper bounds on Sq for square sources and a 5-function source in 2D 
computed from (11). The upper bound for any finite-size source is asymptotically ~ Pe, 

Pe/lnPe in the large Pe limit. In 3D, the distinction 

the 



but for the 5-function source it is 



between cubic sources and a 5- function source is more apparent as shown in Fig 3(b) 
upper bound for a 5-function source behaves ^ 



Pe in 3D. We stress that these mixing 
enhancement bounds apply for any statistically homogeneous and isotropic flows stirring 
sources with these shapes. 



Simulation results for the random sine flow, shown in Fig. 4(a) for 2D and Fig. 4(b) for 
3D qualitatively confirm the behavior of the enhancements suggested by the upper limits. 
As the source size shrinks, the measured mixing enhancement gets smaller in a way that is 
remarkably similar to the bounds. In these simulations Pe is varied by decreasing k at a 
fixed values of L, U and T. Other values of T and other (shorter) wavelengths of the stirring 
flow were also checked, producing similar plots. These 2D simulation results have recently 
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FIG. 3: The theoretical upper bounds for (a) square sources with sizes a = -L/2, -L/10, L/50, and a 
(5-function source; (b) cubic sources with sizes a = -L/5, L/50, -L/500, and a (^-function source (from 
top to bottom). 
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FIG. 4: Measured mixing enhancements for (a) square sources with sizes o = L/2, L/IO, L/50 and 
a 5-function source, (b) cubic sources with sizes a = L/5, L/50, L/500 and a 5-function source 
(from top to bottom). 



been confirmed quantitatively by a PDE computation |32j . 

The simulations also show that the upper estimates can give the correct quantitative 
behavior of as a function of Pe. Indeed, in Fig [5] we plot the upper bound on £^o for the 
5-function source in 3D and the measured enhancement from the simulations. The upper 
bound, which scales anomalously ~ v'Pe, is an excellent predictor of the data. From this 
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FIG. 5: Mixing enhancement for a (5-function source. The soUd hne is the upper bound for any 
SHI flow and the data are mixing enhancements for the random sine flow measured in the discrete 
particle simulations. 



we conclude that the random sine flow is an "almost-optimal" mixer (among statistically 
homogeneous and isotropic flows) for this source-sink distribution. 

V. SUMMARY AND CONCLUSIONS 



We have devised an accurate and computationally eflicient particle method to study 
hydrodynamic variance suppression by a mixing flow. Rigorous upper bounds for the mixing 
enhancement Sq, i.e., the effective diffusion enhancement factor, were compared to measured 
enhancements for the simple random sine flow. A key prediction of the analysis in Refs. [211 
23] is that the source-sink shape is a determining factor in the mixing enhancement of any 
flow. The simulation results reported here show that the upper estimates give the correct 
qualitative picture as regards the Pe and source-shape dependence of Sq. 

Future work should focus on investigating enhancements of other stirring flows. No 
attempt has been made here to flnd a more efficient stirring flow (or, indeed, the most 
efficient flow, if there is one) whose enhancement approaches more closely (or perhaps even 
saturates) the upper bound. It is remarkable that the simple random sine flow appears to 
saturate the upper bound scahng So ~ \/Pe in 3D. 

Stirring with appropriate turbulent solutions to the Navier-Stokes equation is also of 
signiflcant interest. The central question here is, is statistically homogeneous and isotropic 
turbulence generically an efficient mixer? The answer may depend on the source-sink dis- 
tribution. 
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